This script analyzes and plots data for Symbiontic Integration 2021 respirometry data for Pacuta larvae.

Setup

Set up workspace, set options, and load required packages.

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse') 
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car') 
if ("lme4" %in% rownames(installed.packages()) == 'FALSE') install.packages('lme4') 
if ("lmerTest" %in% rownames(installed.packages()) == 'FALSE') install.packages('lmerTest') 
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales') 
if ("cowplot" %in% rownames(installed.packages()) == 'FALSE') install.packages('cowplot') 
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2') 
if ("effects" %in% rownames(installed.packages()) == 'FALSE') install.packages('effects') 
if ("emmeans" %in% rownames(installed.packages()) == 'FALSE') install.packages('emmeans') 
if ("multcomp" %in% rownames(installed.packages()) == 'FALSE') install.packages('multcomp') 

#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('lme4')
library('lmerTest')
library('scales')
library('cowplot')
library('effects')
library('emmeans')
library('multcomp')

Data manipulation

Load data from LoLinR.

PRdata<-read.csv("Pacu2021/Output/Respiration/PR/oxygen_P_R_calc.csv") #load data

Format data columns.

#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]

#format columns
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Plate<-as.factor(PRdata$Plate)
PRdata$Run<-as.factor(PRdata$Run)

Calculate a P:R ratio using gross photosynthesis and respiration calculated as GP:R.

PRdata$ratio<-abs(PRdata$GP.nmol.org.min)/abs(PRdata$R.nmol.org.min) #calculate ratio with absolute values

View outliers. Remove outliers for P:R data.

boxplot(PRdata$R.nmol.org.min)

boxplot(PRdata$P.nmol.org.min)

boxplot(PRdata$GP.nmol.org.min)

boxplot(PRdata$ratio)

PRdata<-PRdata%>%filter(ratio < 4) #filter out outlier for PR ratio
boxplot(PRdata$ratio)

Plotting

Respiration

Plot data as a box plot.

r_plot<-PRdata %>%
    group_by(Temperature)%>%
    ggplot(., aes(x = as.factor(Temperature), y = R.nmol.org.min, group=Temperature)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.2, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); r_plot 

PRdata$Date<-as.factor(PRdata$Date)

r_plot1<-PRdata %>%
    group_by(Temperature, as.factor(Date))%>%
    ggplot(., aes(x = as.factor(Temperature), y = R.nmol.org.min, group=Temperature)) +
    facet_wrap(~Date)+
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.2, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); r_plot1

Plot data as a mean dot plot.

r_plot_b<-PRdata %>%
    group_by(Temperature)%>%
    summarise(mean=mean(R.nmol.org.min), sd=sd(R.nmol.org.min), N=length(R.nmol.org.min), se=sd/sqrt(N))%>%
    ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) + 
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.1, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); r_plot_b

Photosynthesis (Net)

For photosynthesis filter out runs 7 and 8, these runs were outliers that didn’t work due to equipment issues.

Plot data as a box plot.

p_plot<-PRdata %>%
    group_by(Temperature)%>%
    ggplot(., aes(x = as.factor(Temperature), y = P.nmol.org.min, group=Temperature)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); p_plot

p_plot1<-PRdata %>%
    group_by(Temperature)%>%
    ggplot(., aes(x = as.factor(Temperature), y = P.nmol.org.min, group=Temperature)) +
    facet_wrap(~Date)+
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); p_plot1

Plot data as a mean dot plot.

p_plot_b<-PRdata %>%
    group_by(Temperature)%>%
    summarise(mean=mean(P.nmol.org.min), sd=sd(P.nmol.org.min), N=length(P.nmol.org.min), se=sd/sqrt(N))%>%
    ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) + 
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.005, 0.15), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); p_plot_b

Photosynthesis (Gross)

Plot data as a box plot.

gp_plot<-PRdata %>%
    group_by(Temperature)%>%
    ggplot(., aes(x = as.factor(Temperature), y = GP.nmol.org.min, group=Temperature)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); gp_plot

gp_plot1<-PRdata %>%
    group_by(Temperature)%>%
    ggplot(., aes(x = as.factor(Temperature), y = GP.nmol.org.min, group=Temperature)) +
    facet_wrap(~Date)+
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); gp_plot1

Plot data as a mean dot plot.

gp_plot_b<-PRdata %>%
    group_by(Temperature)%>%
    summarise(mean=mean(GP.nmol.org.min), sd=sd(GP.nmol.org.min), N=length(GP.nmol.org.min), se=sd/sqrt(N))%>%
    ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
    geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) + 
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
    scale_y_continuous(limits=c(-0.005, 0.15), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="none",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); gp_plot_b

P:R Ratio

Plot data as a box plot.

pr_plot<-PRdata %>%
    group_by(Temperature)%>%
    ggplot(., aes(x = as.factor(Temperature), y = ratio, group=Temperature)) +
    geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("P:R")))) +
    scale_y_continuous(limits=c(0, 4), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); pr_plot

pr_plot1<-PRdata %>%
    group_by(Temperature)%>%
    ggplot(., aes(x = as.factor(Temperature), y = ratio, group=Temperature)) +
    facet_wrap(~Date)+
    geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
    geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) + 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("P:R")))) +
    scale_y_continuous(limits=c(0, 4), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); pr_plot1

Plot data as a mean dot plot.

pr_plot_b<-PRdata %>%
    group_by(Temperature)%>%
    summarise(mean=mean(ratio), sd=sd(ratio), N=length(ratio), se=sd/sqrt(N))%>%
    ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
    geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
    geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) + 
    geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+ 
    xlab("Temperature") + 
    scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
    scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
    ylab(expression(bold(paste("P:R")))) +
    scale_y_continuous(limits=c(0, 2), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
    theme_classic() + 
    theme(
      legend.position="right",
      axis.title=element_text(face="bold", size=16),
      axis.text=element_text(size=12, color="black"), 
      legend.title=element_text(face="bold", size=14), 
      legend.text=element_text(size=12)
      ); pr_plot_b

Generate final figure

Box plot assembled figure.

full_fig<-plot_grid(r_plot, p_plot, gp_plot, pr_plot, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")

ggsave(filename="Pacu2021/Figures/Respiration/PR/respirometry_fig_boxplot.png", plot=full_fig, dpi=500, width=18, height=6, units="in")

Mean plot assembled figure.

full_fig_b<-plot_grid(r_plot_b, p_plot_b, gp_plot_b, pr_plot_b, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")

ggsave(filename="Pacu2021/Figures/Respiration/PR/respirometry_fig_mean.png", plot=full_fig_b, dpi=500, width=18, height=6, units="in")

Analysis

Respiration

Change data columns to factors.

PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Date<-as.factor(PRdata$Date)
PRdata$Run<-as.factor(PRdata$Run)

View respiration data.

hist(PRdata$R.nmol.org.min)

Build model and examine for respiration across temperatures.

Rmodel1<-aov(R.nmol.org.min~Temperature, data=PRdata) 
summary(Rmodel1)
##              Df  Sum Sq  Mean Sq F value Pr(>F)  
## Temperature   2 0.00617 0.003087   2.537  0.082 .
## Residuals   180 0.21903 0.001217                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check for effects of day.

Rmodel1b<-aov(R.nmol.org.min~Temperature*Date, data=PRdata) 
summary(Rmodel1b)
##                   Df  Sum Sq  Mean Sq F value Pr(>F)  
## Temperature        2 0.00617 0.003087   2.562  0.080 .
## Date               2 0.00454 0.002269   1.884  0.155  
## Temperature:Date   1 0.00123 0.001228   1.019  0.314  
## Residuals        177 0.21326 0.001205                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance.

qqPlot(residuals(Rmodel1))

## [1] 131 123
leveneTest(residuals(Rmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  2.0868 0.1271
##       180

Conduct post hoc test.

emm<-emmeans(Rmodel1, ~Temperature)
pairs(emm)
##  contrast estimate      SE  df t.ratio p.value
##  27 - 31   0.01214 0.00590 180   2.056  0.1022
##  27 - 34  -0.00127 0.00688 180  -0.184  0.9815
##  31 - 34  -0.01340 0.00743 180  -1.805  0.1708
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Net Photosynthesis

Build model and examine for Photosynthesis across temperatures.

Pmodel1<-aov(P.nmol.org.min~Temperature, data=PRdata) 
summary(Pmodel1)
##              Df Sum Sq  Mean Sq F value   Pr(>F)    
## Temperature   2 0.0157 0.007851   7.986 0.000475 ***
## Residuals   180 0.1769 0.000983                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check for effect of day.

Pmodel1b<-aov(P.nmol.org.min~Temperature*Date, data=PRdata) 
summary(Pmodel1b)
##                   Df  Sum Sq  Mean Sq F value   Pr(>F)    
## Temperature        2 0.01570 0.007851   8.328 0.000349 ***
## Date               2 0.01008 0.005040   5.347 0.005565 ** 
## Temperature:Date   1 0.00001 0.000009   0.009 0.922918    
## Residuals        177 0.16685 0.000943                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance.

qqPlot(residuals(Pmodel1))

## [1] 88 58
leveneTest(residuals(Pmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   2  7.2644 0.0009246 ***
##       180                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conduct post hoc test.

emm<-emmeans(Pmodel1, ~Temperature)
pairs(emm)
##  contrast estimate      SE  df t.ratio p.value
##  27 - 31   0.00579 0.00531 180   1.090  0.5215
##  27 - 34   0.02466 0.00618 180   3.989  0.0003
##  31 - 34   0.01888 0.00667 180   2.828  0.0144
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Gross Photosynthesis

Build model and examine for gross photosynthesis across temperatures.

GPmodel1<-aov(GP.nmol.org.min~Temperature, data=PRdata)
summary(GPmodel1)
##              Df Sum Sq  Mean Sq F value Pr(>F)  
## Temperature   2 0.0247 0.012359   4.241 0.0159 *
## Residuals   180 0.5246 0.002914                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check for effect of day.

GPmodel1b<-aov(GP.nmol.org.min~Temperature*Date, data=PRdata)
summary(GPmodel1b)
##                   Df Sum Sq  Mean Sq F value Pr(>F)  
## Temperature        2 0.0247 0.012359   4.395 0.0137 *
## Date               2 0.0253 0.012660   4.501 0.0124 *
## Temperature:Date   1 0.0014 0.001445   0.514 0.4745  
## Residuals        177 0.4978 0.002812                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance.

qqPlot(residuals(GPmodel1))

## [1]  41 131
leveneTest(residuals(GPmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  1.2835 0.2796
##       180

Conduct post hoc test.

emm<-emmeans(GPmodel1, ~Temperature)
pairs(emm)
##  contrast estimate      SE  df t.ratio p.value
##  27 - 31  -0.00635 0.00914 180  -0.695  0.7666
##  27 - 34   0.02593 0.01065 180   2.436  0.0417
##  31 - 34   0.03228 0.01149 180   2.809  0.0152
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

P:R

Build model and examine for PR ratio across temperatures.

PRmodel1<-aov(ratio~Temperature, data=PRdata) 
summary(PRmodel1)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## Temperature   2   2.89  1.4471   4.087 0.0184 *
## Residuals   180  63.73  0.3541                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check for effect of day.

PRmodel1b<-aov(ratio~Temperature*Date, data=PRdata) 
summary(PRmodel1b)
##                   Df Sum Sq Mean Sq F value  Pr(>F)   
## Temperature        2   2.89  1.4471   4.339 0.01447 * 
## Date               2   4.62  2.3078   6.920 0.00128 **
## Temperature:Date   1   0.08  0.0837   0.251 0.61701   
## Residuals        177  59.03  0.3335                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of model for residual normality and variance.

qqPlot(residuals(PRmodel1))

## [1] 159  94
leveneTest(residuals(PRmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   2  4.1833 0.01676 *
##       180                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conduct post hoc test.

emm<-emmeans(PRmodel1, ~Temperature)
pairs(emm)
##  contrast estimate    SE  df t.ratio p.value
##  27 - 31     0.189 0.101 180   1.873  0.1495
##  27 - 34     0.312 0.117 180   2.660  0.0230
##  31 - 34     0.124 0.127 180   0.975  0.5936
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Generate summary

Generate summary of all respiration data.

summary<-PRdata%>%
  group_by(Temperature)%>%
  filter(R.nmol.org.min<0)%>%
  filter(P.nmol.org.min>0)%>%
  filter(GP.nmol.org.min>0)%>%
  summarise(N=length(R.nmol.org.min),
            Mean_R=mean(R.nmol.org.min), 
            SE_R=sd(R.nmol.org.min)/sqrt(length(R.nmol.org.min)),
            Mean_P=mean(P.nmol.org.min), 
            SE_P=sd(P.nmol.org.min)/sqrt(length(P.nmol.org.min)),
            Mean_GP=mean(GP.nmol.org.min), 
            SE_GP=sd(GP.nmol.org.min)/sqrt(length(GP.nmol.org.min)),
            Mean_PR=mean(ratio), 
            SE_PR=sd(ratio)/sqrt(length(ratio)))

summary
## # A tibble: 3 × 10
##   Temperature     N  Mean_R    SE_R Mean_P    SE_P Mean_GP   SE_GP Mean_PR
##   <fct>       <int>   <dbl>   <dbl>  <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 27             81 -0.0728 0.00345 0.0577 0.00358   0.130 0.00531    1.91
## 2 31             54 -0.0828 0.00374 0.0485 0.00367   0.131 0.00580    1.66
## 3 34             32 -0.0738 0.00690 0.0308 0.00295   0.105 0.00799    1.55
## # … with 1 more variable: SE_PR <dbl>
summary%>%
  write_csv(., "Pacu2021/Output/Respiration/PR/mean_respiration.csv")